

******************************************************************
* experiment 2 *
******************************************************************

********************* Table 3
clear
use "data/both.dta", clear
keep if experiment==2
set seed 20251202


* Count sample sizes for each treatment
count if treatment==1
local n_treat1 = r(N)
count if treatment==2
local n_treat2 = r(N)

gen roll = .
replace roll=2 if _n<9
replace roll=1 if _n<5
gen screens =1
replace screens = 2 if mod(_n,2)	
replace screens=. if _n>8
gen mean=.
gen tstat=.
gen ksstat=.
gen pvalue1=.
gen qvalue1 =.
gen pvalue2=.
gen qvalue2 =.
gen pfamily2 = .
gen df=.
gen ci95l =.
gen ci95h =.


ttest roll1==3.47 if treatment==1 // test of mean, JDB
replace pfamily2=r(p) if _n==1
replace mean=r(mu_1) if _n==1
replace tstat=r(t) if _n==1
replace df 		= r(df_t) 			if _n==1
local   df 		= r(df_t) 			
local critical_value = invt(`df', 0.975) // to calaculate the 95% CI that NHB require, but stata doesn't easily report. 
replace ci95h  	=r(mu_1)+ `critical_value'*  r(se) if _n==1
replace ci95l  	=r(mu_1)- `critical_value'*  r(se) if _n==1
ttest roll1==3.47 if treatment==2
replace pfamily2=r(p) if _n==2
replace mean=r(mu_1) if _n==2
replace tstat=r(t) if _n==2
replace df 		= r(df_t) 			if _n==2
local   df 		= r(df_t) 			
local critical_value = invt(`df', 0.975) // to calaculate the 95% CI that NHB require, but stata doesn't easily report. 
replace ci95h  	=r(mu_1)+ `critical_value'*  r(se) if _n==2
replace ci95l  	=r(mu_1)- `critical_value'*  r(se) if _n==2
mgof roll1 = (2/36)*roll1 + 1/36  if treatment==1, mc ksmirnov reps(100000)  // distribution, JDB
replace ksstat=r(ksmirnov) if _n==1
replace pfamily2=r(p_ksmirnov) if _n==3
mgof roll1 = (2/36)*roll1 + 1/36  if treatment==2, mc ksmirnov reps(100000)  // distribution, JDB
replace ksstat=r(ksmirnov) if _n==2
replace pfamily2=r(p_ksmirnov) if _n==4
qqvalue pfamily2, method(hochberg )  qvalue(qfamily2) 
list pfamily2 qfamily2 if qfamily2!=.

replace pvalue1= pfamily2 if _n<3
replace qvalue1= qfamily2  if _n<3
replace pvalue2= pfamily2[_n+2] if _n==1 | _n==2
replace qvalue2=qfamily2[_n+2] if _n==1 | _n==2

    **** FHB, roll 1 ****
gen pfamily3 = .	
ttest roll1==2.53 if treatment==1 // test of mean, JDB
replace pfamily3=r(p) if _n==1
replace mean=r(mu_1) if _n==3
replace tstat=r(t) if _n==3
replace df 		= r(df_t) 			if _n==3
local   df 		= r(df_t) 			
local critical_value = invt(`df', 0.975) 
replace ci95h  	=r(mu_1)+ `critical_value'*  r(se) if _n==3
replace ci95l  	=r(mu_1)- `critical_value'*  r(se) if _n==3


ttest roll1==2.53 if treatment==2 // test of mean, JDB
replace pfamily3=r(p) if _n==2
replace mean=r(mu_1) if _n==4
replace tstat=r(t) if _n==4
replace df 		= r(df_t) 			if _n==4
local   df 		= r(df_t) 			
local critical_value = invt(`df', 0.975) 
replace ci95h  	=r(mu_1)+ `critical_value'*  r(se) if _n==4
replace ci95l  	=r(mu_1)- `critical_value'*  r(se) if _n==4

mgof roll1 = 1/6 if treatment==1, mc ksmirnov  reps(100000) // distribution, JDB
replace pfamily3=r(p_ksmirnov) if _n==3
replace ksstat=r(ksmirnov) if _n==3
mgof roll1 = 1/6 if treatment==2, mc ksmirnov reps(100000)  // distribution, JDB
replace pfamily3=r(p_ksmirnov) if _n==4
replace ksstat=r(ksmirnov) if _n==4
qqvalue pfamily3, method(hochberg )  qvalue(qfamily3) 
format pfamily3 %9.3f 
format qfamily3 %9.3f 
list pfamily3 qfamily3 if qfamily3!=.

replace pvalue1= pfamily3[_n-2] if _n==3 | _n==4
replace qvalue1= qfamily3[_n-2] if _n==3 | _n==4
replace pvalue2= pfamily3 		if _n==3 | _n==4
replace qvalue2= qfamily3		if _n==3 | _n==4

    **** JDB, roll 2 ****
gen pfamily4 = .	
ttest roll2==1.53  if treatment==1 // test of mean, JDB
replace mean=r(mu_1) if _n==5
replace tstat=r(t) if _n==5
replace pfamily4=r(p) if _n==1
replace df 		= r(df_t) 			if _n==5
local   df 		= r(df_t) 			
local critical_value = invt(`df', 0.975)
replace ci95h  	=r(mu_1)+ `critical_value'*  r(se) if _n==5
replace ci95l  	=r(mu_1)- `critical_value'*  r(se) if _n==5

ttest roll2==1.53  if treatment==2 // test of mean, JDB
replace mean=r(mu_1) if _n==6
replace tstat=r(t) if _n==6
replace pfamily4=r(p) if _n==2
replace df 		= r(df_t) 			if _n==6
local   df 		= r(df_t) 			
local critical_value = invt(`df', 0.975) 
replace ci95h  	=r(mu_1)+ `critical_value'*  r(se) if _n==6
replace ci95l  	=r(mu_1)- `critical_value'*  r(se) if _n==6
mgof  roll2 = (2/36)*(5-roll2) + 1/36  if treatment==1 , mc ksmirnov reps(100000)  // distribution, JDB
replace pfamily4=r(p_ksmirnov) if _n==3
replace ksstat=r(ksmirnov) if _n==5
mgof  roll2 = (2/36)*(5-roll2) + 1/36  if treatment==2 , mc ksmirnov reps(100000) // distribution, JDB
replace pfamily4=r(p_ksmirnov) if _n==4
replace ksstat=r(ksmirnov) if _n==6
qqvalue pfamily4, method(hochberg )  qvalue(qfamily4) 
format pfamily4 %9.3f 
format qfamily4 %9.3f 
list pfamily4 qfamily4 if qfamily4!=.

replace pvalue1= pfamily4[_n-4] if _n==5 | _n==6
replace qvalue1= qfamily4[_n-4] if _n==5 | _n==6
replace pvalue2= pfamily4[_n-2] if _n==5 | _n==6
replace qvalue2= qfamily4[_n-2] if _n==5 | _n==6

    **** FHB, roll 2 ****
gen pfamily5 = .	
ttest roll2==2.5 if treatment==1 // test of mean, FHB
replace pfamily5=r(p) if _n==1
replace mean=r(mu_1) if _n==7
replace tstat=r(t) if _n==7
replace df 		= r(df_t) 			if _n==7
local   df 		= r(df_t) 			
local critical_value = invt(`df', 0.975) 
replace ci95h  	=r(mu_1)+ `critical_value'*  r(se) if _n==7
replace ci95l  	=r(mu_1)- `critical_value'*  r(se) if _n==7


ttest roll2==2.5 if treatment==2 // test of mean, FHB
replace mean=r(mu_1) if _n==8
replace tstat=r(t) if _n==8
replace pfamily5=r(p) if _n==2
replace df 		= r(df_t) 			if _n==8
local   df 		= r(df_t) 			
local critical_value = invt(`df', 0.975) 
replace ci95h  	=r(mu_1)+ `critical_value'*  r(se) if _n==8
replace ci95l  	=r(mu_1)- `critical_value'*  r(se) if _n==8


mgof  roll2 = 1/6 if treatment==1 , mc ksmirnov reps(100000) // test of distribution, FHB
replace pfamily5=r(p_ksmirnov) if _n==3
replace ksstat=r(ksmirnov) if _n==7
mgof  roll2 = 1/6 if treatment==2 , mc ksmirnov reps(100000) // test of distribution, FHB
replace pfamily5=r(p_ksmirnov) if _n==4
replace ksstat=r(ksmirnov) if _n==8
qqvalue pfamily5, method(hochberg )  qvalue(qfamily5) 
format pfamily5 %9.3f 
format qfamily5 %9.3f 
list pfamily5 qfamily5 if qfamily5!=.

replace pvalue1= pfamily5[_n-6] if _n==7 | _n==8
replace qvalue1= qfamily5[_n-6] if _n==7 | _n==8
replace pvalue2= pfamily5[_n-4] if _n==7 | _n==8
replace qvalue2= qfamily5[_n-4] if _n==7 | _n==8

order  roll screens mean ci95h ci95l tstat pvalue1 qvalue1 ksstat pvalue2 qvalue2 
keep   roll screens mean ci95h ci95l tstat pvalue1 qvalue1 ksstat pvalue2 qvalue2 
format mean ci95h ci95l tstat pvalue1 qvalue1 ksstat pvalue2 qvalue2  %9.3f
drop if _n>8
texsave  using "output/table3.tex" , replace  


local fn "Note: Sample sizes: Treatment 1 N=`n_treat1', Treatment 2 N=`n_treat2'."
texsave using "output/table3.tex", autonumber varlabels hlines(-2) nofix replace  footnote("`fn'") frag

